/* benchmark guru interfaces */
#include <Utilities/RandomGenerator.h>
#include <Configuration.h>
#include <mpi/communicator.h>
#include <mpi/collectives.h>
#include <Utilities/Timer.h>
#include <OhmmsPETE/OhmmsVector.h>
#include <Message/OpenMP.h>
#include <Numerics/OhmmsBlas.h>
#include <benchmark/codelets.hpp>
#define USE_ALLOCATOR
#if defined(USE_ALLOCATOR)
inline void* allocate(size_t PS, size_t n)
{
  void *d;
  posix_memalign(&d,PS,n);
  return d;
}
#endif

namespace dns_ops
{
  template<typename T, typename RT>
  INLINE0 static void 
    apply_blas(T* restrict Q,  T* restrict Q_tmp, T* restrict Q_old
        , const T* restrict rhs, RT one_third, RT delta, int n)
    {
      BLAS::axpy(n,one_third,rhs,Q);
      BLAS::copy(n,Q_old,Q_tmp);
      BLAS::axpy(n,delta,rhs,Q_tmp);
    }
}

using namespace APPNAMESPACE;

int main(int argc, char** argv)
{
  mpi::environment env(argc,argv);
  mpi::communicator mycomm;
  OhmmsInfo ohmms("rk4",mycomm.rank(),0,1);

  typedef  double real_type;
  const int n_var=3;
  int niters=10;
  int nz=12288;
  int nc=576;
  bool debug=false;

  size_t mypagesize=sysconf(_SC_PAGESIZE);
  int ic=0;
  while(ic<argc)
  {
    string a(argv[ic]);
    if(a.find("-d")<a.size()) //get the z
      nz=atoi(argv[++ic]);
    if(a.find("-n")<a.size()) //number of pencils
      nc=atoi(argv[++ic]);
    else if(a.find("-i")<a.size()) 
      niters=atoi(argv[++ic]);
    else if(a.find("-p")<a.size())
      mypagesize=atoi(argv[++ic]);
    else if(a.find("-debug")<a.size())  
      debug=true;
    ++ic;
  }

  size_t from_system=sysconf(_SC_PAGESIZE);

  //mypagesize=std::min(from_system,mypagesize);
  cout << "system pagesize = " << sysconf(_SC_PAGESIZE) <<endl;
#if defined(USE_ALLOCATOR)
  cout << "allocator pagesize = " << mypagesize <<endl;
#else
  cout << "using PETE Vector. " << endl;
#endif

  double t0=0,t1=0,t2=0,t3=0;
//  int np=omp_get_max_threads();
//#pragma omp parallel for reduction(+:t0,t1,t2,t3)
//  for(int ip=0; ip<np; ++ip)
#pragma omp parallel reduction(+:t0,t1,t2,t3)
  {
    int ntot=n_var*nz*nc;
#if defined(USE_ALLOCATOR)
    real_type *Q = static_cast<real_type*>(allocate(mypagesize,ntot*sizeof(real_type)));
    real_type *rhs = static_cast<real_type*>(allocate(mypagesize,ntot*sizeof(real_type)));
    real_type *Q_tmp = static_cast<real_type*>(allocate(mypagesize,ntot*sizeof(real_type)));
    real_type *Q_old = static_cast<real_type*>(allocate(mypagesize,ntot*sizeof(real_type)));
#else
    Vector<real_type> Q(ntot);
    Vector<real_type> rhs(ntot);
    Vector<real_type> Q_tmp(ntot);
    Vector<real_type> Q_old(ntot);
#endif

    double one_third=1.0/3.0;
    double delta=0.01/2.0;
    Timer myclock;
    vector<double> timing(4,0.0);
    for(int i=0; i<niters*4; ++i)
    {
#if defined(USE_ALLOCATOR)
      myclock.restart();
      dns_ops::apply(Q,Q_tmp,Q_old,rhs,one_third,delta,ntot);
      timing[0]+=myclock.elapsed();

      myclock.restart();
      dns_ops::apply1(Q,rhs,one_third,ntot);
      dns_ops::apply2(Q_tmp,Q_old,rhs,delta,ntot);
      timing[1]+=myclock.elapsed();

      myclock.restart();
      dns_ops::apply_blas(Q,Q_tmp,Q_old,rhs,one_third,delta,ntot);
      //do_rk4(Q,Q_tmp,Q_old,rhs,&one_third,&delta,&ntot);
      timing[2]+=myclock.elapsed();
#else
      myclock.restart();
      dns_ops::apply(Q.data(),Q_tmp.data(),Q_old.data(),rhs.data(),one_third,delta,ntot);
      timing[0]+=myclock.elapsed();

      myclock.restart();
      dns_ops::apply1(Q.data(),rhs.data(),one_third,ntot);
      dns_ops::apply2(Q_tmp.data(),Q_old.data(),rhs.data(),delta,ntot);
      timing[1]+=myclock.elapsed();

      myclock.restart();
      dns_ops::apply_blas(Q.data(),Q_tmp.data(),Q_old.data(),rhs.data(),one_third,delta,ntot);
      timing[2]+=myclock.elapsed();

      //using expression template
      myclock.restart();
      Q += one_third*rhs;
      Q_tmp=Q_old+delta*rhs;
      timing[3]+=myclock.elapsed();
#endif
    }

#if defined(USE_ALLOCATOR)
  free(Q);
  free(rhs);
  free(Q_tmp);
  free(Q_old);
#endif

    t0=timing[0];
    t1=timing[1];
    t2=timing[2];
    t3=timing[3];
  }

  {
    vector<double> timing(4,0.0);
    timing[0]=t0; timing[1]=t1; timing[2]=t2; timing[3]=t3;

    mpi::reduce(mycomm,timing);

    {//write out for confluence wiki
      const double rk4=1.0;
      double nodefac= rk4/static_cast<double>(mycomm.size()*niters*omp_get_max_threads());
      ostringstream os;
      os << "MPI Nodes = " << mycomm.size() 
        << " OpenMP Threads = " << omp_get_max_threads() << endl;
      os << "fft-dim num-ffts mpi omp ops loop1 loop2 f77 pete niters mpi omp" << endl;
      os << nz << " " << nc << " " 
        << mycomm.size() << " " << omp_get_max_threads()  << " "
        << nz*nc*n_var;
      for(int i=0; i<timing.size(); ++i) os << "  " << timing[i]*nodefac;
      app_log() << os.str() << " " << niters  << " "
        << " RK4 " << endl;
    }
  }
  return 0;
}

